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We consider the cross-correlation search for periodic gravitational waves and its potential appli¬ 
cation to the low-mass x-ray binary Sco X-1. This method coherently combines data not only from 
different detectors at the same time, but also data taken at different times from the same or different 
detectors. By adjusting the maximum allowed time offset between a pair of data segments to be 
coherently combined, one can tune the method to trade off sensitivity and computing costs. In 
particular, the detectable signal amplitude scales as the inverse fourth root of this coherence time. 

The improvement in amplitude sensitivity for a search with a maximum time offset of one hour, 
compared with a directed stochastic background search with 0.25-Hz-wide bins is about a factor 
of 5.4. We show that a search of one year of data from the Advanced LIGO and Advanced Virgo 
detectors with a coherence time of one hour would be able to detect gravitational waves from Sco 
X-1 at the level predicted by torque balance over a range of signal frequencies from 30 to 300 Hz; if 
the coherence time could be increased to ten hours, the range would be 20 to 500 Hz. In addition, 
we consider several technical aspects of the cross-correlation method: We quantify the effects of 
spectral leakage and show that nearly rectangular windows still lead to the most sensitive search. 

We produce an explicit parameter-space metric for the cross-correlation search, in general, and as 
applied to a neutron star in a circular binary system. We consider the effects of using a signal tem¬ 
plate averaged over unknown amplitude parameters: The quantity to which the search is sensitive 
is a given function of the intrinsic signal amplitude and the inclination of the neutron-star rotation 
axis to the line of sight, and the peak of the expected detection statistic is systematically offset 
from the true signal parameters. Finally, we describe the potential loss of signal-to-noise ratio due 
to unmodeled effects such as signal phase acceleration within the Fourier transform time scale and 
gradual evolution of the spin frequency. 


I. INTRODUCTION 

The low-mass x-ray binary (LMXB) Scorpius X-1 (Sco 
X- m is one of the most promising potential sources 
of gravitational waves (GWs) which may be observed 
by the generation of GW detectors—such as Advanced 
LIGOP], Advanced Virgo[3] and KAGRA[3] —which will 
begin operation in 2015 with the first Advanced LIGO 
observing run, and Advanced Virgo and KAGRA obser¬ 
vations expected to follow in the coming years. Sco X-1 
is presumed to be a binary consisting of a neutron star 
which is accreting matter from a low-mass companion; its 
parameters are summarized in Table |lj Nonaxisymmetric 
deformations in the neutron star can give rise to gravi¬ 
tational radiation, most of which is emitted at twice the 
rotation frequency of the neutron star[Tn]0 Such defor¬ 
mations can be maintained by the accretion of matter 
onto the neutron star. It has been conjectured m that 
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^ Additionally, unstable rotational modes of the neutron star, or r 
modes m can lead to GW at 4/3 of the neutron star’s rotational 
frequency. 


the neutron star’s rotation may be in an approximate 
equilibrium state, where the spin-up torque due to ac¬ 
cretion is balanced by the spin-down due to gravitational 
waves. Scorpius X-l’s high x-ray flux implies a high ac¬ 
cretion rate, which makes it the most promising potential 
source of observable GWs among known LMXBsfT^. 

Since Sco X-1 is not seen as a pulsar, its rotation fre¬ 
quency is unknown. There is also residual uncertainty 
in the orbital parameters which determine the Doppler 
modulation of the signal, monochromatic in the neutron 
star’s rest frame, which reaches the solar-system barycen- 
ter (SSB). This parameter uncertainty limits the effec¬ 
tiveness of the usual coherent search for periodic grav¬ 
itational waves [TO]. The first search for GW from Sco 
X-1 with the first generation of interferometric GW de¬ 
tectors, using data from the second LIGO science run[8], 
was limited to six hours of data for this reason. A subse¬ 
quent search with data from the fourth LIGO science run 
m used a variant of the cross-correlation method devel¬ 
oped to search for stochastic GW backgrounds, treating 
Sco X-I as a random unpolarized monochromatic source 
with a known sky location[TS]j^ 


^ Other methods have been developed, specialized to search for 
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TABLE I. Parameters of the low-mass x-ray binary Scorpius 
X-1. Since the sky position is determined to microarcsecond 
or better accuracy, the relevant astrophysical parameters with 
residual uncertainty are those describing the orbit. Those are 
the projected semimajor axis ap = a sin i of the neutron star’s 
orbit, the orbital period Porb, and the time tasc at which the 
neutron star crosses the ascending node (moving away from 
the observer), measured in the solar-system barycenter. The 
orbital eccentricity of Sco X-1 is believed to be small[T], and 
the present work presumes the orbit to be circular for sim¬ 
plicity; consideration of eccentric orbits would add two search 
parameters which are determined by the eccentricity and the 
argument of periapseO |6]. Note that the observational con¬ 
straint in [1] is not on Up itself, but on the radial velocity 
amplitude Ki = of the primary. We could have formu- 

•^orb 

lated the parameter space in terms of and Porb rather than 
ap and Porb, but this has no significant impact on the accu¬ 
racy of the method, since the uncertainty in ap is dominated 
by that associated with K\. Finally note that the orbital ref¬ 
erence time tasc (which we quote as the time of ascension of 
the compact object, 1/4 cycle before the time of inferior con¬ 
junction of the companion quoted in 0 ) can be propagated 
to a later epoch by adding an integer number of periods, at 
the cost of increasing the uncertainty due to the uncertainty 
in the period itself. 


Parameter 

Value 

Reference(s) 

Right ascension 

16‘’19“55.0850" 

[S] from 

Declination 

-15°38'24.9" 

[8] from [9] 

Distance (kpc) 

2.8 ±0.3 

m 

Up (sec) 

1.44 ±0.18 

[S] from [T] 

Gsc (GPS sec) 

897753994 ± 100 

m 

Porb (sec) 

68023.70 ± 0.04 

m 


The stochastic analysis formed the inspiration for a 
new method to search for periodic gravitational waves 
with a model-based cross-correlation statistic which takes 
into account the signal model for continuous GW emis¬ 
sion from a rotating neutron star[5T]. (This method has 
also been adapted |22] to search for young neutron stars 
in supernova remnants.) The present work further devel¬ 
ops some of the details of this method and the specifics 
of applying it to search for gravitational waves from Sco 
X-1 and, by extension, other LMXBs. 

The paper is organized as follows: Section [IT] reviews 
the basics of the method and the construction of the com¬ 
bined cross-correlation statistic using a new, streamlined 
formalism. Section m works out the statistical proper¬ 
ties of the cross-correlation statistic, including the first 
careful determination of the effects of signal leakage and 


LMXBs. These include summing over contributions from side¬ 
bands created by Doppler modulation llfilITTl . searching for such 
modulation patterns in doubly-Fourier-transformed data ll8lll9l . 
and fitting a polynomial expansion in the Doppler-modulated 
GW phase [20l. 


the unknown value of the inclination angle of the neutron 
star’s axis to the line of sight. It also considers in detail 
how the sensitivity of the model-based cross-correlation 
search should compare to the directed unmodeled cross¬ 
correlation search for a monochromatic stochastic back¬ 
ground. Section |IV] considers two effects related to the 
dependence of the statistic on phase-evolution param¬ 
eters such as frequency and binary orbital parameters: 
a systematic offset of the maximum in parameter space 
from the true signal parameters (which depends on the 
unknown inclination angle), and the quadratic falloff of 
the signal away from its maximum. The latter is en¬ 
coded in a parameter space metric, which we construct 
in general as well as for the LMXB search both in its 
exact form and in limiting form relevant if the obser¬ 
vation time is long compared to the orbital period. In 
Sec.0we consider limitations to the method from inaccu¬ 
racies in the signal model, either due to slight variations 
in frequency (“spin wandering”) arising from an inexact 
torque-balance equilibrium, or due to phase acceleration 
during a stretch of data to be Fourier transformed. Fi¬ 
nally, in Sec. |VI| we summarize our results and consider 
the expected sensitivity of this search to Sco X-1. 


II. CROSS-CORRELATION METHOD 


The cross-correlation method is derived and described 
in detail in |21j . In this section, we review the fundamen¬ 
tals, using a more streamlined formalism and including a 
more careful treatment of signal-leakage issues and nui¬ 
sance parameters. 


A. Short-time Fourier transforms 


Because the signal of interest is nearly monochromatic, 
with slowly varying signal parameters, it is convenient to 
describe the analysis in the frequency domain by dividing 
the available data into segments of length Tgft and calcu¬ 
lating a short-time Fourier transform (SFT) from each. 
Since the sampling time St is typically much less than 
the SFT duration Tgft, we can approximate the discrete 
Fourier transform of the data by a finite-time continuous 
Fourier transform. If we use the index K to label both the 
choice of detector and the selected time interval, which 
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has midpoint tx, the SFT will 


N-l 

XKk = ^ xxitx — 'Fsft/2 + j St) 

j=0 

+Tsft/2 

J iic —T'sft/2 

^iic+Tsft/2 

= (- 1 )'=/ , 

( 2 . 1 ) 

where the frequency corresponding to the fcth bin of the 
SFT is 


fk = kSf = 


Tsft 


( 2 . 2 ) 


In practice, the data are often multiplied by a window 
function wj = w ) before being Fourier trans¬ 

formed, so that (|2.ip becomes 


N-l 

j=0 

rtK+Tstt/2 /f_f 
(- 1 )'= / 


Tsft 


In this work we assume that the windowing function is 
nearly rectangular with some small transition at the be¬ 
ginning and end, so that leakage of undesirable spectral 
features is suppressed, but the effects of windowing on 
the signal and noise can otherwise be ignored. The im¬ 
plications of other window choices are considered in Ap¬ 
pendix 1^ 


B. Mean and variance of Fourier components 


ani 


C0 


E[nK{t)nL{t')]=SKL r d/ . 

( 2 . 6 ) 

If we write the noise contribution to the SFT labeled by 
K as 

N-l 

nxk = riKj St 

3=0 ^2.7) 

ftK+Tsft /2 

« / n;f(f)e"‘^’"b-[iK-Trft/2])/*, 

tjf —'Att/2 


then (2.51 implies E[nKk] = 0 and we can use (2.6) to 
show that 


E [n_R•/c?^z,t] ~ ^KL Ski Tsft 


Sxifk) 


( 2 . 8 ) 


(As detailed in Appendix]^ this is not the case for non¬ 
trivial windowing, where noise contributions from differ¬ 
ent frequency bins are correlated.) If we can estimate 
the noise PSD Sxifk), we can “normalize” the data to 
define (as in [53]) 


(2.3) 


ZKk — XKk 


TsftSx 


which has mean 


E [ZKk] = t‘'Kk = hxk 1/ rri c 

V -tsftOK 

unit covariance 

E [{ZKk — fJ-Kk){ZLl — fJ-Ll)*] = SkL Ski 
and zero “pseudocovariance” 

E [(zKk — tkKk){zLl — IJ-Ll)] = 0 


(2.9) 


(2.10) 


( 2 . 11 ) 


( 2 . 12 ) 


Let the data 


(This is because the real and imaginary parts of each zxk 
are independent and identically distributed.) 


XK{t) = hxit) + nxit) 


(2.4) 


in SFT K consist of the signal hx (t) plus random instru¬ 
mental noise nx (t) with one-sided power spectral density 
(PSD) Sx{\f\) so that its expectation value is 


E[nx{t)] = 0 


(2.5) 


C. Signal contribution to SFT 

The signal from a rotating deformed neutron star is 
determined by various parameters of the system, which 
can be divided into the following categories uni. 


^ Note that the factor e appears in Eq. (2.25) of [2j with 

the wrong sign in the exponent. However, given for integer 
k, this phase correction is simply the sign (—1)^ so the complex 
conjugate does not change it. 


Strictly speaking, we should allow for data from adjacent 
SFT intervals in the same detector to be correlated, but 
we assume that the autocorrelation function Knit — t') = 
f^oo df falls off quickly compared to Tgft, so 

that we can neglect the correlation between noise in different 
time intervals. 
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(i) Amplitude parameters: intrinsic signal amplitude 
ho, the angles l and '0 which define the orientation 
of the neutron star’s rotation axis (t is the incli¬ 
nation to the line of sight and tp is a, polarization 
angle from celestial west to the projection of the 
rotation axis onto the plane of the sky), and the 
signal phase <i>o at some reference time. 

(ii) Phase-evolution parameters: intrinsic phase evolu¬ 
tion (frequency and frequency derivatives) of the 
signal, as well as parameters such as sky location 
and binary orbital parameters which govern the 
Doppler modulation of the signal. 

Those parameters determine the signal received by a 
gravitational-wave detector at time t as 

h{t) = ho {F^A+ cos$(<) -I- FxAx sin4>(t)) (2.13) 

where id|_ and Fx are the antenna pattern functions [TOl 
[21] which change slowly with time as the Earth rotates. 
The signal contribution to a SET can be estimated by 

hif (t) « ho{F^A+ cos($if -b 2tt fxit - Ik]) 

+ FxAx sin($if -b 27r/if [t - tie])} 

where we have Taylor expanded the phase about the time 
Ik'- 


-\-2TTfK{t — tx) ■ (2.15) 


The validity of this approximation will be one of the lim¬ 
iting factors which determines the choice of SET duration 


Tgft, as detailed in Sec. VB 


The form of (2.141 includes the following parameters 
and definitions: 


(i) A+ = '' and Ax = cos i depend on the incli¬ 

nation L of the rotation axis to the line of sight. 

(ii) The antenna patterns F^ and F^ depend on the 
detector in question, the sidereal time at tx, the 
sky position a, 6, and the polarization angle ip. 

(hi) The relationship t{t) between the SSB time and the 
time at the detector depends on the sky position 
and timej^ Thus the phase 4’(t(t)) depends on time, 
detector, 4>o, /o,/i,..., sky position and—in the 
case of a binary—the binary orbital parameters. 


The signal contribution to bin k of SET K is 


hxk 


A , — \ A 

ho(-l)'=e‘^^ + + ^ ^ ' 


■<^Trft(/fe - /if) 

(2.16) 


® Specifically, if fhet is the position of the detector and k is the unit 
vector pointing from the source to the SSB, t(f) ^ t — f\jet ' k/c. 



_I__^__I__^__I_ 

-2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 


(/ - f')/Sf 


FIG. 1. Plot of if—f') defined in (2.171 which determines 
the signal contribution to a given frequency bin of a short 
Fourier transform (SFT) of duration Tstt according to (2.161. 
Since the spacing between frequency bins is Sf — 1/Tsft, there 
will be, for a given signal frequency fx, one bin whose value 
of KKk = ifk — fK)/Sf lies between each pair of vertical solid 
lines. 


where we have defined 


ptK+Tsit/2 

ifk -fK)= / g-i2.(/.-/rr)(t-trr) 

—T’sft/2 

= Tsft sinc([/fc - /ifjTsft) 

(2.17) 


in terms of the normalized sine function sine a = 

7ra 

This is plotted in Fig. ^ The signal contribution will 
be largest in the fc^fth Fourier bin, defined by 


kx := 


/if 

Jf 


L/if?"sftl 


(2.18) 


whose frequency /j,^ is closest to fx- (We have intro¬ 
duced the notation that [a] is the closest integer to a.) 
It will prove useful to define, similarly to [23] 

Kxk = k — /ifTsft = = kx + (k — kx) (2.19) 

6f 

where 

kx = =kK- /ifTsft , (2.20) 

of 

so that — 5 < kx < 5 . A simple search would consider, 
from each SFT K, only the Fourier component 


® Previous sensitivity estimates |21l I22| noted that (0) = ^sft 
and therefore replaced each of the finite-time delta functions with 
the SFT length Tgft, but a more careful treatment requires that 
we keep track of spectral leakage caused by the signal frequency 
not being centered in a SFT bin. 

^ Note that our definition of Kxk differs by a sign from the one 
used in m- 
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closest in frequency to the signal frequency Jk at the 
search parameters. However, as we will see, the sensi¬ 
tivity of the search can be improved by including contri¬ 
butions from additional adjacent bins, so we indicate by 
K-k the set of bins to be considered from SFT K, and 
we will construct a detection statistic using XKk for all 

k e JCk- 

We can then write 


hKk « smc{KKk)e''^‘^ 


which means that, from (2.101 



Tsit 

( 2 . 21 ) 


E [zKk] = fkKk 

rpK A _ • j^K A 

« ho(-l)'= sinc(KKk)E‘^^ + ^ ^ ^ 


2Tsft 


Sk 

( 2 . 22 ) 


D. Construction of the cross-correlation statistic 


For a given choice of signal parameters, which deter¬ 
mine kk for each SFT, and therefore kkr for each Fourier 
component, it is useful to deling 


ZK = 


(-1)^ sinc(KKfc)zj^fc 
= ^ (-l)''sinc(K_R-fc)z_R-fc 


(2.23) 


k^Kh 


This is still normalized so that 

E [{zk — ^kK){zL — ^J-L)*] = 5 kl 
E [{zK - iji-k){zl - Ml)] = 0 


where now 


Mic ~ hoe 


A , A 

i$ic V ^ 


X „ / 2rsft 

K\ 


Sk 


(2.24a) 

(2.24b) 


(2.25) 


If we define vectors indexed by SFT number, we can write 
(2.24) and (2.25) in matrix form as 


E 


f; [z] = M 

(z-m)(z-m)^ 


= 1 


E (z — m)(z — =0 


(2.26a) 

(2.26b) 

(2.26c) 


Note that computations can be made more efficient by use of the 
identity sinc(Kiffc) = so (-1)'= smc(Kxfe) = 

(_l)fejc sin(7rKjy) where only the final factor depends on the 
bin index k S Afx- 


where 1 is the identity matrix, 0 is a matrix of zeros, 
(•)'^' indicates the matrix transpose and (•)^ the matrix 
adjoint (complex conjugate of the transpose). 

A real cross-correlation statistic p can be constructed 
by defining a Hermitian matrix W and constructing p = 
z^Wz = Tr(Wzzl). [Our chosen form of W will be 
defined in (2.35).] Equation (2.26) tells us that 


E [zzl] =1-1- M/i,! (2.27) 

where the second term is a matrix with elements 

2Tsft 


2" = 


FkFl = h^'E^K'^Le 


"Fl-l 




(2.28) 


where /S.^kl = — 4>l is the difference between the 

modeled signal phases in the two SFTs and F^-l is a 
geometrical factor which depends on l and ip as follows 
[compare Eq. (3.10) of pT]]: 

Tkl = I {F^E^Al + E^E^Al 

+ i[E^E^-F^E^]A+A^) 

= 1 ( -4^ (^iC^L ^ 

4 V 2 


-f iA+Ax (a^b^-b^a^) 


uKL-, 


[{a^a^-b^b^) cos4^ 


+ {a^b^ + h^a^) sin4V'] 


(2.29) 


where we have used the fact that the tp dependence of 
the antenna patterns Fp^ can be written in terms of 
the amplitude modulation (AM) coefficients and b^ 
as 


F^ = cos2ilj + b^ sin2ilj (2.30) 

F^ = — sin 2?/; -|-6^ cos 2^ . (2-31) 

The AM coefficients [To] are determined by the relevant 
sky position, detector and sidereal time. They can be 
defined pS] as and b^ = where e“^ 

and are a polarization basis defined using one basis 
vector pointing west along a line of constant declination 
and one pointing north along a line of constant right as¬ 
cension. Note that t and ip are properties of the source 
which do not change for different SFT pairs, while 
and depend only on the SFT (detector and sidereal 
time) and sky position. It is also useful to note that the 
combinations 

F^Fip + Ff 4 = + b^b^ = 10F|J'2 (2.32a) 

F^F^ - F^Fip = a^b^ - b^a^ = lOF^i; (2.32b) 

are independent of ip. 
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Since terms in change signs if we vary cos i and ip, 
which are unknown, it is convenient, as proposed in |21j . 
to work with the average over those quantities, which 
picks out the “robust” part: 

^Tl = + b^b^) (2.33) 

Note that is real and non-negative, while Tkl is 
complex. On the other hand, Tkl can be factored into 
IkIl^ while cannot. If we define (again as in 

[23| . but with a different overall normalization) “noise- 
weighted AM coefficients” and by dividing by 
\l construct Tkl from those, we can write 


^^K^J'L = = H^Ckl 


(2.34) 


or, as a matrix equation, = h^G. Note that [2T] 
did not consider issues of spectral leakage responsible for 
S/f, and used a different convention for the placement of 
complex conjugates in atomic cross-correlation term, so 
their Qkl would be equal to in the present nota- 

Gl 


tion. Similarly, our 


from 


corresponds to the combination 


V'’ 

As noted in m, an “optimal” combination of cross¬ 
correlation terms would use a weight W proportional 
to G. However, as described above, we work with 
Gkl = in order to avoid specifying 

the parameters cos i and ip. For reasons of computa¬ 
tional cost to be detailed later, we limit the possible 
set of SFT pairs KL included in the cross-correlation 
to some set V, in particular by requiring that K < L 
and — til < Tmax- Then we define the Hermitian 
weighting matrix W by 

'NG^jPI KLeV 

Wkl = { Nid’^^D* LKer (2.35) 

0 otherwise 

so that the cross-correlation statistic is 
p = z^Wz = Tr(Wzz^) 

= ^ E {gTl^k^l + GTl^zkzI) 

= Ar E "^KL E E (-l)^~^sinc(Kjcfc)sinc(KL^) 

KLgV keICK t&K.L 


X {e^^'^^^ZKkZLi + e 




ZKkZLi) 


(2.36) 


Since we assume that the list of pairs V includes no 
autocorrelations, the matrix W contains no diagonal el¬ 
ements which implies Tr(W) = 0. We will later in¬ 
troduce, and use when convenient, the notation that a 
labels a (nonordered) pair of SFTs KL S V. 


III. STATISTICS AND SENSITIVITY 

In this section we consider in detail the statistical 
properties of the cross-correlation statistic p which were 
sketched in a basic form in |21j . In particular, we con¬ 
sider the impact on the expected sensitivity of spectral 
leakage and unknown amplitude parameters, and com¬ 
pare the sensitivity of a cross-correlation search to the 
directed stochastic search by analogy to which it was de¬ 
fined. 


A. Mean and variance of cross-correlation statistic 


The expectation value of the cross-correlation statistic 
is 


E[p]=E [Tr(Wzzt)] = Tr(W) + hi Tr(WG) 

=/igTr(WG) =/x^W /2 

where we have used the fact that W is traceless. The 
variance is 

Var(p) = E [p2] _ E [pf = E [z^WzzWz] - (/x^W/x)^ 

(3.2) 

The first term can be evaluated by writing z = (z—/x)-t-/x; 
after some simplification we have 


Var(p) =E 
+ 


(z — /x)^W(z — 

2/x^WV 


/x)(z-/x)^W(z-/x) 


(3.3) 


Ordinarily we would need to know something about the 
fourth moment of the noise distribution to evaluate the 
expectation value, but since W contains no diagonal ele¬ 
ments, and the different elements of z—/x are independent 
of each other, the expectation value can be evaluated us¬ 
ing only the variance-covariance matrix of z to give 

Var(p) = Tr -p 2/x^WV = Tr -p 2hl Tr W^G 

(3.4) 

We choose the normalization constant N so that p has 
unit variance in the limit hi —>■ 0, i.e., 

l = Tr(W2)=^^W^fiWLx = 2iV2 ^ 

K L KL^V 

(3.5) 


® Note that eq (3.10) of [2l]is also missing a factor of (—1)^^ 
which should appear in ^LkE' omission was pointed 

out in |22| . but eq (5) of |22| included the wrong sign in the phase 
correction and failed to stress that the relevant frequency is 
rather than /jy- 


Note that if we analogously constructed the matrix to include 
only diagonal terms, i.e., constructed a statistic only out of auto¬ 
correlations, the statistic would be equivalent to that used in the 
PowerFlux method|26|. 
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i.e., 


N-^ = 2J2 \gTl =2 5] 


klgv 


"2 ■z:2 


klgv 


{^TlY (3.6) 


Written in terms of SFT pairs, the expectation value 
of the statistic is 


E[p] = hlTTiWG) 

= Nhl Y. {gtAl + GTlGkl) 

KLeV 

= Nhl2 Y S^5iric™RerKL 


(3.7) 


KLeV 


Looking at (2.291 we see that the real part oiVKL has a 
piece proportional to and a piece that depends on 
ijj-. 


ReFifi = - 


5Al 


■Al 


-pave 

KL' 


A‘^ — A‘^ 


(3.8) 

The sum over SFT pairs KL can be broken down as a 
sum over detector pairs, over time offsets Ik — tL, and 
over the time stamp ^(t/c + ti) halfway between the 
time stamps of the SFTs in the pair. In an idealized 
long observing run, if the detector noise is uncorrelated 
with sidereal time, the sum over + tL) means we 

are averaging the two expressions {a^a^ + and 

{a^a^ + b^b^)iFfF^ - F^F^) (the latter of which de¬ 
pends on the polarization angle ip) over sidereal time. 
Because the former is positive definite and the latter is 
not, this average tends to suppress the V'-dependent term. 


_l_^ 

This is in addition to the fact that — 


> 


possibly substantially, depending on the value of t, as il ¬ 
lustr ated in Fig. If we neglect the second term in (3.81, 
(13.71) becomes 


F[p\^Nhl- 


bAl 


—-2 


E 

KL£V 


’='2 -=>2 


( f^ave \ 

iiCLj 


= K)^ 



(3.9) 


where 


Ff = ho 


5 Al + Al 


(3.10) 


is the combination of ho and cos t that we can estimate 
by filtering with the averaged template. 

Since we have normalized the statistic so that Var(/9) = 
1 for weak signals, the expectation value ( |3.9[ ) is an ex¬ 
pected signal-to-noise ratio for a signal with a given 
This means that if we define a SNR threshold such 
that p > p^'^ corresponds to a detection, the signal will 
be detectable if 


hf>A 



- 1/4 


(3.11) 



FIG. 2. Plot of AlfPjl the coefficients of the 


Al+At 


two contributions to RePxL in (3.8 1 . The factor 
is also equal to | where is the combination of 

ho and cos t approximately measured by the cross-correlation 
statistic, as shown in, e.g., Eq. (3.91 


B. Impact of spectral leakage on estimated 
sensitivity 

Finally, we consider the impact of the leakage factors 
of the form sinc^(Kiyfc) on the expectation 

value. Expanding these expressions, we have 


F[p]^{hf)^i 2 Y (rKl) 

^ KL^V 

X Y, sinc^(Kiyfc) Y, sinc^(KLf) 






1/2 


(3.12) 


If we choose only the “best bin” kx = kx from each 
SFT, defined by (2.181), we have 


^2 


^ = sinc^(K*:) 


(3.13) 


If, instead of the best bin whose frequency is closest 
to /if, we take the m closest bins to define JCx, the sum 
becomes 


"2 

“if 


fejc+L(m-l)/2j 

= Y, sinc^(Kiffe) = Y, sinc^(«;if/c) 

k&K-K fe=fcjf-|'(m-l)/2] 

[(m-l)/2j 

= 'Y^ sinc^(kif-h s) 

s=-\(m-l)/ 2 \ 

(3.14) 


where [aj < a and \a\ > a are the integers below and 
above a, respectively. Note that, because of the iden- 





































TABLE II. Contributions to (H^), defined in (3.181, from in¬ 
clusion of multiple SET bins. We see that using a single bin 
from each SET leads to only around 77.4% of the maximum 
sensitivity given by (3.151, but that we can recover over 90% 


of this sensitivity by using two bins and over 95% by using 
four bins from each SET. This table applies for rectangularly 
windowed data; using other window options further reduces 
the expected SNR, as described in Appendix The table 
also assumes that the various Doppler modulations move the 
signal frequency around to accomplish an average over the 
fractional offset of the signal frequency from the center of the 
bin. The validity of this approximation is explored in m- 


m 

1 

2 

3 

4 

5 

6 

Contribution 

0.774 

0.129 

0.028 

0.019 

0.009 

0.007 

Cumulative 

0.774 

0.903 

0.931 

0.950 

0.959 

0.966 


tit3p]X]^-oo sinc^(K -I- s) = 1, valid for any k, the best 
we can do by including more bins is < 1 and there- 
forJ^ 


E[p]<ihff 2 Y. (r-)^ 


(3.15) 


klgv 


The sensitivity associated with the inclusion of a finite 
number of bins from each SET will depend on the value of 
— 5 < kk < 5 corresponding to the signal frequency /k 
in each SFT. We can get an estimate of this by assuming 
that, over the course of the analysis, the Doppler shift 
evenly samples the range of k values, and writing 


E [p] « {hff j2{E?y Y 

V KL^V 


= {hfr{E^)i2 


E 


(3.16) 


KLeV 


with 


^) = 


sinc^(K + s) 


(3.17) 


where (•)« indicates an average over the possible offsets 
within the bin. We can numerically evaluate 




= ^(sinc"(^c + s))^ 

S 

[(m-l)/2j 1/2 

= Y, ^ / sinc^(K-b s) dK (3.18) 

pmj2 

= 2 / sinc^AcdK 
Jo 


as shown in Table m 

Since most cross-correlation searches will be compu¬ 
tationally limited, the question of how many bins to in¬ 
clude from each SFT is one of optimization of resources. 
The value of E [p] for a given /iq®, and therefore the 
sensitivity of the search, can be increased by including 
more frequency bins from each SFT, but this will involve 
more computations and therefore more computational re¬ 
sources. If instead those resources were put into a search 

with a larger Tmax, the value of Y1,kl^v would 

be higher. Naively, one might expect the computing cost 
to scale with the number of terms to be combined, and 
therefore with the square of the number of bins taken 
from each SFT. So increasing from m = 1 to m = 2 
could take up to 4 times the computing cost. On the 
other hand, for a fixed number of bins, we suppose that 
the cost will scale with the number of SFT pairs to be 
included times the number of parameter space points to 
be searched. Typical behavior will be for the density of 
points in parameter space to scale with some 

integer value of d; as described in Sec. |IVB[ for a search 
over frequency and two orbital parameters of an LMXB, 
as long as T^ax is small compared to the binary orbital 
period, d = 3. Since the number of SFT pairs at fixed 
observation time will also scale like Tmax, the overall com¬ 
puting cost will scale like , and quadrupling the com¬ 
puting time would mean multiplying the possible 


and thus the number of terms in the sum (3.161 by . 
This would increase E [p\ for a given dg® by a factor of 
42 ( 3 + 1 ) = 23+T. For d = 3, this is 2^/'* « 1.19, which 
is very slightly more than the benefit « 1.17 from 
including a second bin from each SFT. However, the as¬ 
sumption that computing cost scales like w? is likely an 
overestimate (since most of the operations can be done 
once per SFT rather than once per pair), so it is generally 
advisable to use at least two bins from each SFT. 


11 This 

. 1/2 


1 - 1 / 


writing 

E oo 

s= —o 


sinc(«: + s) = 

pi2'Ks(t — t') _ 


is most easily proved by 
^2 dt and using 

E“_oo< 5 (i-i' + s). 

1^ Previous sensitivity estimates |21II22| were missing the factor of 
and therefore slightly overestimated the sensitivity. 


C. Sensitivity estimate for unknown amplitude 
parameters 


The cross-correlation statistic is 
Var(p) « 1 and, according to (3.16) 


normalized so that 
and now adopt¬ 


ing the notation that a refers to an unordered allowed 
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pair of SFTs, 


E[p] = ihr) 


effN2/"2 



= (3.19) 


where /iq® is the combination of hg and cos t given in 
(3.10), and is a property of the search which can 


be determined from noise spectra, AM coefficients, and 
choices of SFT pairs, without knowledge of signal param¬ 
eters other than the approximate frequency and orbital 
parameters. Even if the noise in each data stream is 
Gaussian distributed, the statistic, which combines the 
data quadratically, will not be. It was observed in m 
that each individual cross correlation between SFTs is 
Bessel distributed; the of the optimal sum is consid¬ 
ered in Appendix both in its exact form and a nu¬ 
merical approximation. For simplicity, in what follows 
we assume that the central limit theorem allows us to 
treat the statistic as approximately Gaussian, with mean 
and unit variancef^ 

We consider the sensitivity estimates in [^, which 
implicitly assume the values of t and '0 are known and 
used to construct the expected cross correlation used in 
weighting the terms in the statistic. [In our notation 
this wo uld m ean using Gkl rather than in the def¬ 


inition (2.35) of W.] Here we perform the analogous 


calculation, assuming we’re using in the construc¬ 
tion of the statistic. Thus the probability of exceeding a 
threshold p*® will be 

poo 

P(p > = / f{p\ho,i-,tp)dp 


exp [p- dp 

1 _ 1 - hlg{c) 


1 

^20 


= - erfc 
2 


72 


= - erfc 
2 


72 


(3.20) 


where 


5 Ai + Ai „„„ 5 


P(0 ~ —^7™ = f^(l + 6cos^ i -b cos"^ 

(3.21) 

The threshold associated with a false alarm probability 
a is 


= 72 erfc ^(2a) 


(3.22) 


but the sensitivity /iq®"® associated with a false dismissal 
probability /? will now be defined, following a procedure 


analogous to the one in [29] , by marginalizing over the 
unknown inclination i (since we have neglected the 0 
dependence in E MS 


i-/3 = p(p>7®iho = hr^) 

= (p(p>7®l/^o = /^^^bV’)) 


cos 


^ 1 / / 7"-(fer^)^g(0 

2\ V 72 


(3.23) 


So to get a se nsitiv ity estimate, we need to find the /iq®' 
which solves (3.23), i.e., 


2(1-0) 


erfc ^ — 


pth (^^6118)2 ^ave 5 


(1 -I- 6cos^ 6 -I- COS^ i) 


72 72 16 

= erfc ^erfc" 72 a) - >5®®^ [l + 6x^ -b dx 


so that the approximate sensitivity is 


= 


/5®®72 


= (5' 


■eff\-2/.-2\2 


i:(f"0 


(3.24) 


- 1/4 


(3.25) 

Equation (3.24) defines 5®® as a specific function of a 


and 0, so the approximate sensitivity correction due to 
marginalizing over cos t can be worked out independently 
of the details of the search. We show some sample values 
Table for a and 0 values between 1% and 10%, and 
also for single-template a values corresponding to overall 
false alarm probabilities in the same range, assuming a 
trials factor of 10®. We see that the ho sensitivity is 
modified by between 39% and 67% in these cases. 


D. Scaling and comparison to directed stochastic 
search 


We consider here the behavior of ( 3.25[ ) (or equiva¬ 
lently ( 3.19[ )) with parameters such as the observing time 
Tobs and allowed lag time Tmax, which is effectively a co¬ 
herence time. As noted in EH, the detectable ( |3.25[ ) 
scales like one over the fourth root of the number of SFT 


Note that this approximation is less accurate in the tails of the 
distribution. Unfortunately, for a search over many independent 
templates, the most interesting statistic will necessarily be in the 
tails. For example, with 10® templates, even a 1% false alarm 
probability for the loudest statistic value would correspond to 
a single-template false alarm probability of 10“^®. See|28| for 
specific examples of this. 


Note that if we had kept the i/j-dependent term in l |3.8[ >, the 
resulting E [p] /Hq would depend not only on both l and but 
also on the detector geometry and pairs of SFTs and a numerical 
solution to the equivalent of i 3.23| l would have to be performed 
anew for basically each sensitivity estimate. 
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TABLE III. Approximate modification of search sensitivity, as a fnnction of desired false alarm probability a (corresponding 
to a statistic threshold of p**') and false dismissal probability /3, resnlting from filtering with a template averaged over the 
signal parameters cos t and ip. (The second set of a values is chosen to correspond to interesting single-template false alarm 
probabilities with a trials factor of 10®.) The detectable signal amplitude ho'’"® (3.251 is proportional to V5®*^. The table shows, 
for a variety of choices of a and /3, how the corrected factor \/5®® calculated according to (3.241 compares to the standard 
expression S = erfc“^(2a) -|- erfc“^(2/3) which would apply from filtering with known values of the parameters cost and ip. 
Note that using the worst-case value cos t = 0 shows that 1 < S°^/S <3.2. 


a 



5 

P 



P 



P 


0.10 

0.05 

0.01 

0.10 

0.05 

0.01 

0.10 

0.05 

0.01 

0.10 

1.3 

1.81 

2.07 

2.55 

3.49 

4.45 

6.27 

1.39 

1.47 

1.57 

0.05 

1.6 

2.07 

2.33 

2.81 

4.15 

5.16 

7.03 

1.42 

1.49 

1.58 

0.01 

2.3 

2.55 

2.81 

3.29 

5.42 

6.52 

8.47 

1.46 

1.52 

1.60 

10“® 

6.0 

5.15 

5.40 

5.89 

12.73 

14.16 

16.40 

1.57 

1.62 

1.67 

5 X 10“^° 

6.1 

5.23 

5.48 

5.96 

12.96 

14.40 

16.64 

1.57 

1.62 

1.67 

10-10 

6.4 

5.40 

5.66 

6.14 

13.48 

14.93 

17.20 

1.58 

1.62 

1.67 


pairs included in the sum Ejg 


Lsens 
flQ 


= (5^S)-2(s2)2iVp,i,,((rr)^ 


- 1/4 




2 /coff\-2/"2\2 


4(r 


KLJ 


SkSl 


- 1/4 


(3.26) 


search. If we only allow simultaneous pairs of SFTs, the 
number of pairs included in the sum (3.25) becomes 


^simul 

'pairs 


Ndet{Ndet — 


1 ) 


Fobs 

Fsft 


(3.29) 


which makes the signal strength to which the search is 
sensitive 


The approximate number of pairs for a search of data 
from A^det detectors, each with observing time Fobs (so 
that the total observation time is iVdetFobs)) with maxi¬ 
mum lag time F^ax > T^ft is 


N, 


pairs 


-'''det 


Fobs Fjoax 

Fsft Fgft 


so the sensitivity scaling is 


(3.27) 


- (^iVd\tFobsFa,ax(5^®)-"(s2)^ (^iSr)) 

(3.28) 

We wish to compare this sensitivity to that of the di¬ 
rected stochastic search (also known as the “radiometer” 
method) defined in [15] and used to set limits on grav¬ 
itational radiation from Sco X-1 [H EO]- The directed 
stochastic search is also an optimally weighted cross¬ 
correlation search, but only includes contributions from 
data taken by different detectors at the same time. We 
first consider the sensitivity of a cross-correlation search 
using our method with this restriction, and then relate 
this to the sensitivity of the actual directed stochastic 


Note that the averages here are not the weighted averages intro¬ 
duced in Sec. |IV[ 


I sens\simul 


(hri 

NdetiNdet - l)FobsFsft(0-'(S2)2 


4(rg'2)- 

SrSl 


- 1/4 


1 - 


1 

Ndet 


Fsft 


- 1/4 


(3.30) 


The directed stochastic search is not quite the same as 
this hypothetical cross-correlation search with simulta¬ 
neous SFTs, however. Most of these differences are ir¬ 
relevant or produce effectively identical calculations. For 
instance, since the Ato, appearing in ( 4.17[ ) is zero for si¬ 
multaneous SFTs, the phase difference A$q, = 27r/oA(i„ 
just encodes the difference in arrival times at the two de¬ 
tectors. Likewise, while the stochastic search assumes a 
random unpolarized signal rather than the periodic signal 
from a neutron star with unknown parameters, this has 
the same effect as our choice to use F|)^£ as the geometri¬ 
cal weighting factor. In fact (as noted in [^) 
is, up to a normalization, the overlap reduction function 
for the directed stochastic search. The one signihcant 
difference is that, since the stochastic search does not 
model the orbital Doppler modulation, it does not have 
access to the signal frequency Jk corresponding to SFT 
K, and therefore cannot localize the expected signal fre¬ 
quency to a bin of width 5f = Thus, instead of 


T,ft ■ 


the optimal combination described by (2.23) or (2.36), it 


must sum with equal weights the contributions ZKkZ^k 
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across a coarse frequency bin of width A/ > -p^fo (see 


(IV B 21 for the definitions of the binary orbital param¬ 
eters relevant to Doppler modulation)!^ The effect is 
to increase the variance of the cross correlation due to 
noise by ^ = A/ Tsft (since there are A/ bins being 
combined, only one of which contains a significant signal 
contribution) so that 


sens\stoch 


(hrn 


AT (AT -.'I'^obs .eoffN-2 /4(r^|,)^ 

A^det(iVdet- 1)^(5 ) 


T sens 
^0 


/’-' 2\-2 


0- 


1 - 


Ndet 


AfT^ 


-1/4 


-1/4 


(3.31) 


The appearance of the factor containing (S^) in the com¬ 
parison is because the directed stochastic search, by com¬ 
bining a larger range of frequency bins, as well as tech¬ 
niques such as overlapping windowed segments, avoids 
some of the usual leakage issues. On the other hand, if 
A/ is chosen to maximize the sensitivity for a given fre¬ 
quency, there will be similar issues with part of the signal 
falling outside the coarse bin at the extremes of Doppler 
modulation. 


To insert concrete numbers, (3.311 tells us that for a 


search with data of equivalent sensitivity from three de¬ 
tectors, a cross-correlation search with T^ax = 3600 s and 
(S^) = 0.9 would provide an improvement in ho sensitiv¬ 
ity over a directed stochastic search with A/ = 0.25 Hz 
of a factor of about 5.4j^ This is consistent with the per¬ 
formance of the two searches in the Sco X-1 Mock Data 
Challenge [35], in which the cross-correlation method was 
able to detect signals with ho almost an order of magni¬ 
tude lower than those detected by the directed stochastic 
method. 


This was not the original motivation for the coarse frequency 
bins in the stochastic cross-correlation pipeline; see for exam¬ 
ple ED, but it has this effect when using the method to search 
for monochromatic signals from neutron stars in binary sys¬ 
tems. Note also that it is sufficient to perform a single sum 
^Kk^Lk across the coarse bin rather than a double sum such 
because, while the frequency bin containing 
the signal is not known, it will be the same bin for both detec¬ 
tors because the unknown phase shift due to the orbit is the same 
for simultaneous SFTs. 

This does not include the fact that the directed stochastic 
method includes a relatively coarse search over frequency, while 
the model-based cross-correlation method must search over many 
more points in frequency and orbital parameter space, as de¬ 
scribed in Sec. |IV ^ This seemingly significant increase in trials 
factor turns out to be swamped by the gain in sensitivity. In 
the comparison above, the same signal will generate a factor of 
almost 30 larger rho value in the cross-correlation search. On the 
other hand, the p threshold to achieve a 5c7 false alarm probabil¬ 
ity would need to be increased only from 5 to 7.8 to overcome a 
trials factor of 10®. Additionally, the search over signal param¬ 
eters in the cross-correlation method allows estimates of those 
parameters. 


Note that, unlike the model-based cross-correlation 
search, the stochastic search is not computationally lim¬ 
ited, with year-long wide-band analyses being achievable 
on a single CPU|35|. Additionally, since it does not as¬ 
sume a signal model (beyond sky localization and ap¬ 
proximate monochromaticity), it is robust against un¬ 
expected features such as orbital parameters outside the 
nominally expected range. However, its sensitivity is fun¬ 
damentally limited by its ignorance of orbital Doppler 
modulation, with a maximum effective coherence time of 


J_ < Parb 
Af ~ 2TTapfo 


(100 Hz A 

V /o ) 


75 sec. 


IV. PARAMETER SPACE BEHAVIOR 


So far we have implicitly ass umed the parameters used 
to construct the signal model (2.16), other than the am¬ 


plitude parameters /iqi cos 6, and i/’, were known when 
constructing the weighted statistic. In order to deter¬ 
mine the phase evolution of the signal, and therefore 
and /ic, we need various phase-evolution parame¬ 
ters {Ai}. (For example, for a neutron star at a known 
sky location with a constant intrinsic signal frequency /o 
in a binary orbit, these are /o and any unknown binary 
orbital parameters.) A slight error in these would lead 
to the appearing in /x and that used to construct W 
being slightly different. In this case we need to go back to 
(3.7) and distinguish between the true A^if l and the one 
assumed in the construction of the filterj^ If we write 
these parameters as {Ai}, let the parameters assumed in 
constructing p be Ai and the true parameters of the sig¬ 
nal be Af. Let A$f^^ and A^kl be the phase difference 
constructed with the true signal parameters and 
the parameters assumed in W, respectively. The effect 
will be t o redu ce the expected SNR E [p\ from the value 
given in (3.19) which it would attain with Ai = Af. The 


modified value is 


E[p] « hlN{^^) 


r. 




r 


* p —' 


(4.1) 


It is also possible for and/or ExSl to differ from their 

assumed values, e.g., if the search parameters include sky po¬ 
sition which can change the amplitude modulation coefficients, 
or a change in Doppler modulation affects the location of the 
signal frequency within the bin. We follow the usual procedure 
of focusing on the dominant effect, which is the change in the 
expected signal phase, and thereby obtain a “phase metric” for 
the cross-correlation search. 
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Now, for Xi close to Af, 

P — A^a) _|_ p* —i(A^»^ — A^'a) 

i ctd ~r i Q,t: 

= 2 Rer„ cos(A$® - A$„) - 2Imra sin(A$* - A$„) 
« 2 Rer„ (^ 1 - i(A$„-A$*) 2 ^+ 2 Imf„(A$„-A$*) 

(4.2) 

if we write the phase difference as 
A$„ - A$« « ^ A$„,,(A, - A|) 


A. Systematic parameter offset 


The result (4.4) not only tells us how the expected SNR 


falls off when the parameters {Ai} used in constructing 
the statistic differ from the true signal parameters {A®}, 
it also shows that the maximum of E [p\ is not actually 
at the signal point Xi = Af, but at the point Xi = Af" 
defined by 


+ 9 ^ A$Q,,ij(Ai - Xi){Xj - Xj) 




(4.3) 


where A$a_i = we obtain, to second order in the 

parameter difference, 

E[p] « (ho)2A(s2) |^2^frRer„^ 

1 - ^ef(A. - Af) -J29.jiX. - Af)(A, - A®) 


where 


2 ^„rrimr„A$„,, 


2 EorrRer„ 

and the parameter space metric is 

1 2 Ea rf,™ (Rc f„A$„,,A$„ j + Im f„A$„,y 


(4.4) 

(4.5) 


9ij = 


2 EarrRer„ 


(4J) 

If we once again neglect the '(/'-dependent piece of Re 
as well as the second derivative term in the metric, we 
have 




1 EkLG'p(“^“^ +b^b^)‘^A^a,i^^a,j 

1 pave A* A(F, 1 

2 2 ^ 


jKLgV ^ a 

where (•)a indicates a weighted average with weighting 
factor [recall fff'l oc ] and 


e,- 


2A+Ax 

“ Al + Al 

EiCLGp(«'^«^ +b^b^)ia^b^ -b^a^)A^KL,^ 


2A+Ai, Ea^r^E"A$„,^ 


Ai + Al 


Ea (rr)^ 


(4.8) 


0 = e® + V2g„(A--A®) 


J 


i.e., at 


\m _ \s _ 

\ \ 2^ ^9ij 


(4.9) 


(4.10) 


where { 3 ^^- } is the matrix inverse of the metric {gij}- 
If the metric is approximately diagonal, so that g~2 ~ 
^, then the offset of the true signal parameters from the 
maximum value of E [p\ is 


1 e 

\ s \ m _ 

Aa — Aa — — — 


2.4+^>, 


E pavepcir 
a Q! ^ a 


A$„ 


2 ga Al + Al 


Ea (rr)' 


A4>„,A4>„ 


(4.11) 

This offset depends on the (generally unknown) value of 
the inclination angle i via A+ = '' and Ax = cos t. 

In particular it has the opposite sign for b € (0, 7 r/ 2 ) and 
b € (7r/2,7r). For a signal detection with unknown i, this 
will have the effect of a systematic error in the measure¬ 
ment of the phase-evolution parameters {Ai}. (Of course, 
one could perform a subsequent analysis which would 
produce an estimate of b, such as a coherent followup of 
the signal candidate, or a cross-correlation search using 
iT^g in place of in the construction of W.) 


B. Parameter space metric 

We return now to consideration of the metric defined 


(4.7) by ( |4.7| ), 


9v = 


lE[p]^ 


2 E[p] 




^(A$„.,A$„.,E ■ (4-12) 


1. Comparison to standard expression for metric 

We can relate this to the usual notation for the phase 
metric. [See, e.g., Eq. (5.13) of |33], which was also used 
inl 22 l .1 




(4.13) 


Note, first of all, that while the standard definition 
of the parameter space metric defines the mismatch as 
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the fractional loss in signal-to-noise squared, our cross¬ 
correlation statistic p is actually the equivalent of what 
is usually called p^. This is because it is quadratic in the 
signal (as is the T statistic, and its expectation value is 
proportional to Hq). 


The connection between (4.12) and (4.13) is made by 


noting that the averages in (4.131 are over data segments, 
while the expression in (4.12) is a weighted average over 
SFT pairs, where the weighting factor is (T^™)^. We can 
relate the two in the special case where the set of pairs 
V contains every combination of SFTs (e.g., by choosing 
Tmax to be the observing time), and by neglecting the 
influence of the weighting factor in the cross-correlation 
metric. In that case, the average can be written as a 
double average over SFTs K and L: 


- 2 


which is just (4.13). Note that this identification can only 


be made in the case where the cross correlation includes 
all pairs of SFTs (or all pairs within some time stretch). 
With a restriction such as \tK — tL\ < one must 

consider the weighted average over pairs, not separate 
averages over SFTs. 


2. Metric for the LMXB search 

We now consider the explicit form of the parameter 
space metric for a neutron star in a circular binary sys¬ 
tem, assuming a constant intrinsic frequency /q. Al¬ 
though the actual values of phase = ‘h(t(<if)) and 


frequency = 


dt 


t^tK 


used via (2.15) to con¬ 


struct the expected cross correlation Gkl include rela¬ 
tivistic corrections, it is sufhcient for the purposes of con¬ 
structing the parameter space metric to limit attention 


= ‘ho + 27r/o tx — 


along the propagation direction from the source. 
(Note that this depends on the detector, but also 
on the time tx-) 


(ii) Qp = 


is the projected semimajor axis of the 


binary orbit, in units of time. 

(iii) Porb is the orbital period of the binary. 

(iv) tasc is a reference time for the orbit, defined as 
the time, measured at the solar-system barycenter, 
when the neutron star is crossing the line of nodes 
moving away from the solar system. 


If we use the identity 


(4.14) 


, 'A + B\ fA-B 

sm A — sin B = 2 cos —-— sin 


we have 


(4.16) 


A$a = 27 r/o<^ Ata - Ad, 


„ . r-Ata 

— 2 Gj) sm —— cos 




orb 


27r - 

{ta ^asc) 


P 


orb 


(4.17) 


where we have defined AdxL = dx — dL, AtxL =tx—tL, 
and txL = ■ 

Note that AdxL will be much less than AtxL unless 
the SFTs K and L are simultaneous. (This is because 
the duration of a SFT will be long compared to the 
light travel time between detectors on the Earth, and 
the Earth’s motion is nonrelativistic.) 

We can now calculate the derivatives appearing in 

(|4T^: 


dA^ 

dfo 


= 27 r< At„ — Adn 


„ . TrAta 

— 2a„ sm —- cos 


27r 

Porb 


(fa ^asc) 


(4.18a) 


dUr, 


A f ■ 

-47r/o sm cos 


27r 

Porb 


(foi ^asc) 


(4.18b) 


= $0 + Stt/o <tx - dx - ttp sin 


where we have defined the following: 


US 


aA$„ 

Stt^/og 

1 ^orb ■ ^ \ 

'-'^asc 

-^orb 


c J 

27T 

p i^K ^asc ) 


dTr/oOf 

sin 

[ dPorb 

-^orb 


sm ■ 


Porb 


sm 


27r 

Porb 


{ia ^asc) 


(4.18c) 


(4.15) 


(i) dx = the projected distance, in seconds, 

from the solar-system barycenter to the detector. 


. irAta . 

X sm- sm 

Porb 

irAta nAta 


27r 

Porb 


(fa ^asc) 


27r - 

(ia ^asc) 


Porb 


(4.18d) 
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3. Approximation for long observation times 


It is relatively simple and straightforward to construct 
the phase metric for a given observation; calculate the 
derivatives (4.18) for each SFT pair and then insert them 
into the weighted average (4.12). However, we can gain 
insight into the behavior of the metric if we consider an 
approximate form which should be valid if the observ¬ 
ing time (e.g., one year) is long compared to the orbital 
period of the LMXB (e.g., 6.8 x 10"*^ s « 19hr for Sco 
X-lpQ[7]). Since the orbital period is not commensurate 
with any of the relevant periods of variation such as the 
sidereal or solar day [the former being relevant for 
and the latter for the noise spectra], it is reasonable to 
assume that — ^asc) samples all phases roughly 

equally, and therefore 


Fa cos 


Fa cos 


= { Fa sin 
27r 


27r ^ ^ 

^asc) 

-^orb 
27r 


-^orb 


Porh 
(j^a ^asc) 


{ia ^asc) 
27T 


sin 


-Porb 


= 0 

(j^a ^asc) 


Fry COS^ 


27r ^ ^ 

^asc) 


= ( Fry sin^ 


27T - 

■ ^asc ) 


orb 


(4.19a) 

= 0 

(4.19b) 

(4.19c) 


— 2 


5tascPorb 


/ {ta)^-t^sc \ / . 2 TTAta \ 

Pit. \ Pori. Pori, /a 


gParbParb = 




p4 

“^orb 


{{ta - tasc)\ (sin^ 


(4.20h) 


irAto 


orb 


-ip 

■^orb 


Atlcos^ 


F 


orb 


(4.20i) 


The metric is not diagonal, but we can neglect the off- 
diagonal elements if 


(gij) ^9ii9jj ■ (4-21) 

One can show that {gfoa^Y < 5/o/o and 

(5/oPo.b)^ < 5 / 0/0 5P„.bPorb as long as 

{{Ata - Ada)\ » al (4.22) 

which should be the case; for Sco X-1, Up = 1.44 s mm- 
Note also that, as long as we include cross correlations 
between non simultaneous SFTs, {{Ata — ^da)‘^)a ~ 
{{Ata)^)^ because the detectors are moving much slower 
than the speed of light. 

We will also have {ga^p^.bf < 9aj,a^ 5P„,bPorb as long as 
the square of the typical time lag Ata is much less than 
{(ta — fasc)^)jj, which will be the case if the maximum 
allowed time lag is much less than the length of the run. 
We can see this by considering the {(ta — ^asc)^)^; if we 
define 


where Fa is any expression not involving ta- 


We then have metric components, from (4.12), of 


5/0/0 = 27 r^ ((A/a - Ada)\ + (sin^ ) 


. 2 f / • 2 7 '’A/q 
5/oap = 47r foap { sm 


5/oP„rb = - 4 


-Porb 

/ , TrAta ttA/o 

^ ( Ata sm - cos ■ 

orb 


Porb Porb 


A 2 r2 / • 2 

9a^a^ = 47r /o ( sm 


orb 


5/otaBC — gapt^sc - 0 

iTT^/oOp / . ttA/q, nAta 

5apP„,b =- m — i sm — cos ■ 

orb 


9iasc ^asc 


foap / . 2 7rA/a 


PL Pori. 


(4.20a) 

(4.20b) 

“(4-20c) 

(4.20d) 

(4.20e) 

I 

“(4-20f) 

(4.20g) 


then 


fJ-T = {ta)c 


L = {{ta - IJ-tLc 


(4.23) 


(4.24) 


should be on the order of the square of the duration of the 
run. In particular, for a run of duration Tobs during which 
the sensitivity of the search remains roughly constant. 


1 


Fobs J- 


fPobs/Z rp2 

edt = . 


Tobs/2 


12 


(4.25) 


But 


{{ta - tasDa = L + (MT “ /asc)^ > L (4-26) 
This leaves only the ratio 


( 5 t, 


L ^ ((^“)a ^osc) 


gtp 


c5PorbPo, 


{{ta /asc)^)^ 


(A^T /asc) 

ctI -I- {g-T - /asc)^ 


(4.27) 

Whether or not this can be neglected seems to come 
down, then, to whether the reference time /asc falls dur¬ 
ing the run. If it falls outside the run, {gp — /asc)^ ^ L 
and the off-diagonal metric element gt^^^p^^b cannot be 













































15 


ignored. However, it is always possible to replace one 
reference time tasc with another = tasc + nPorh sep¬ 
arated by an integer number n of cycles, and thus it is 
always possible to arrange for — t'ascY — Pmh ^ 
and thus obtain an approximately diagonal metric. This 
comes at a cost, though, since there will be a contribu¬ 
tion to the uncertainty in the new reference time due to 
the uncertainty in the orbital period. If the uncertainties 
in the orbital period and the original reference time are 
independent, the uncertainty in the new reference time 
will be given by 

(A4j2 = (At_)^+n2(AP,rb)^ 

= (A4.c)^ + ftasc-^4sc)^ (Ap^^^)2 (4-28) 

orb 

This will become the dominant error if 

l^asc ~ ^asc| > ^orb • (4.29) 

^-^orb 


For Sco X-1, using the parameter uncertainties from [7] 
(see Sec. VI), this is about 


- X 68023.70 s « 5yr 

0.04 ^ 


(4.30) 


Since the tasc quoted in [7] (chosen to minimize their 
Atasc) corresponds to June 2008, this will be the case for 
any GW observations using Advanced LIGO and/or Ad¬ 
vanced Virgo data, unless additional Sco X-1 ephemeris 
updates are made. 

Subject to the aforementioned approximations, the 
metric can be treated as diagonal with non-negligible el¬ 
ements 


the variation of (F^™)^ from pair to pair. Subject to this 
approximation, we have 







ancC!] 



nAta 

Poih 


a 


2T, 


sm' 


9 , 

^ dt 


max J—T^ 


Porh 


1 — sine 


2T 

max 

J^orb 


(4.33) 


(4.34) 


where once again sine a; = ■ Note that this is only 

a rough approximation, since increasing the time offset 
Ata between a pair of SFTs from the same instrument 
(or from well-aligned instruments like the LIGO detec¬ 
tors in Hanford and Livingston) will tend to decrease 
the expected cross correlation as the detectors are ro¬ 
tated out of alignment with each other. We confirm this 
by comparing the approximate expressions to more ac¬ 
curate values calculated using the geometry of the LIGO 
and Virgo detectors and the sky position of Scorpius X-1, 
in Fig. 

Note that some care needs to be taken when compar¬ 
ing our metric expressions to those in [5]. For example. 


combining (4.31a) with ( 4.33[ ) gives us gf^fg « 27r^-“, 
which seems at odds with the analogous expression in 
e.g., Eq. (61) of [B], where the corresponding metric ele¬ 
ment is TT^ . The difference is that the semicoherent 
search in [6] is defined by combining distinct coherent 
segments of length AT, which makes the mean squared 
difference 


9fofo 

g Qip dp 


Staaciaac 

ffPorbPorb 


167r^/(f4 


P^ 


orb 


4"rb 


(4.31a) 

(4.31b) 

(4.31c) 

(4.31d) 


The quantities (At^Q, and ^sin^ j which appear in 

the parameter space metric are constructed by a weighted 
average over SFT pairs. If we consider a search which 
includes all pairs up to a maximum time lag of Tmaxi the 
parameter space resolution, and therefore the required 
number of templates, will depend on Tmax- We can get a 
rough estimate on this dependence by assuming that we 
can write 


{f{Ata))a 


2Trr 


fit) dt 




(4.32) 


which assumes Tobs ^ Tma,x ^ T/ft so that we can replace 
the sum over specific lags with an integral, and neglects 


pAT pAT 

/ / [t — t'Ydtdt' 

Jo Jo 


(AT)2 


I 


fAT pAT-\At\/2 


(AT)2 

rAT 


|At|/2 


(At)^ dt dAt 


1 


(AT)2 


/ lAl 

(At)^(AT- |At|)dAt 

-AT 


= (=- = )(AT)2 = 1(AT)2 (4.35) 


whereas our maximum lag rule — < T^ax gives a 


19 


Note that for Tmax "C -Porb (coherent integration times small 
compared to the binary orbital period), the factor ^sin^ ^ 


tends to (so the number of templates in each direction 

•^-^orb 

grows like the coherent integration time), while for Tmax ^ Porb? 
coherent integration times long compared to the binary orbital 
period, it tends to a constant so the growth in number of tem¬ 
plates in the ap and iasc directions saturates. This is analogous 
to an effect described in m- 
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FIG. 3. Plot of weighted averages and (sin^ ) appearing in the metric components (|4.31|) as a function of 

^ \ / a. _ _ 

maximum allowed lag time Tmax- The dotted lines show the approximate values (4.331 and (4.341 neglecting the variation of 
the weighting factor. The solid line (labeled HLV) shows the value for a search using detectors at the LIGO Hanford, LIGO 
Livingston, and Virgo sites, assuming a source at the sky position of Sco X-1, and that all detectors have the same sensitivity at 
the relevant frequency, and all sidereal times are evenly sampled. The dashed line (HL) shows the same thing for a search using 
only the LIGO detectors at Hanford and Livingston. The actual weighted averages (and therefore the number of templates 
needed to cover the parameter space) are less than the approximate ones, because the geometrical factor (P)),''®)^ weights smaller 
lag times more. 


mean square time difference 


A. Spin wandering 


.Tob, pin(t'+r„„„.TobO _ ,^2 

JO Jmax(i'-Tn,ax.O) ' ' 


rTobs rmin(t'+T„ax.robs) 

Jo Jmax(t'—T„ax,0) ' 

pAT aToba-|Ai|/2, 


I-AtIiaZ~2 

rAT fT>bB-|At|/2 
J-ATJ|At|/2 acaiAl 


(2/3)TobsTLx - mm 


4 

max 


2Tobs'Tinax — T^e 


rj^Z 


(4.36) 


where the assumption Tmax ^ Tobs gives us the result 

(|43^. 


We have assumed so far that the LMXB is in approx¬ 
imate equilibrium, where the spin-up torque due to ac¬ 
cretion is balanced by the spin-down due to gravitational 
waves. Even if this is true on average, the balance will 
not be perfect, and the spin frequency will “wander”. 
This means that rather than a constant frequency fo ap¬ 
pearing in (4.15), there will be a time-varying frequency 

/(t), where t = t — 


is the time measured 


in the neutron star’s rest frame. Thus the phase differ¬ 
ence between SETs K and L will be, rather than just 
= 27r/o[t_R- — ti]. 


V. IMPLICATIONS OF DEVIATION FROM 
SIGNAL MODEL 

So far, we have assumed that the underlying signal 
model contained in (2.21|), along with the phase evolu¬ 


tion (4.15) is correct, although some of the parameters 
may be unknown. We consider two effects which vio¬ 
late this assumption, and their potential impacts on the 
expected SNR (3.19). These are (1) spin wandering, in 


which the frequency is not a constant /o but varies slowly 
and unpredictably with time and (2) the impact of higher 
terms in the Taylor expansion of <i)(i(t)) ab out t = Ik, 
which are neglected in the linear phase model (2.15). The 


former effect will place a potential limit on the coherence 
time Tmax by providing an intrinsic limit to the frequency 
resolution, whereas the latter will constrain our choice of 
SET length T^ft in order that neglected phase accelera¬ 
tion effects not cause too much loss of SNR. 


A^true 


r^L 


= 4>iy - $Z, = 27r 


f{t)di (5.1) 


liK 


We can consider the loss of SNR due to the existence of 
spin wandering, compared to what we would expect if the 
frequency truly were constant. Qualitatively, there are 
two reasons for loss of SNR: first, on short time scales, the 
change in frequency could disrupt the coherence between 
the two SETs in a pair being cross correlated; second, on 
longer time scales, the spin could wander enough that the 
SNR is distributed over different frequency templates. 

To quantify the loss of SNR we follow a calculation 

in (|4T1) and (|T2) 


analogous to that in Sec. IV 
obtain 


e.g.. 


to 


E[pf^^^-E[p] 


E[p] 


ideal 


(5.2) 


where {■)a is a weighted average over SET pairs 
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with weighting factor as before. To estimate 

— A$q,))^^ we assume that the wandering is 

slow enough that we can expand /(t) in a Taylor series 
about Ikl = (tic + tL)/2: 


/(<•) ~ fi^KL) + tifi) 

min(iif, ti) < t < max(iif, ti) (5.3) 

Then 


rti 


Aci>‘^“ - A$kl = 27r 


[/(t) - /o] dt 


'iK 


27r ( [/{Ikl) - M^^kl + /(tifi) 


(Ai^i)' 


(5.4) 


where AIrl = tif — ti Subject to reasonable ass ump tions 
about the randomness of the spin wandering, (5.2) can 
be written in the form 


E[pf^^^-E[p] 


E[p] 


ideal 


^27r2([/(t„)-/o]2)„((At„)2)^ 

+ y([/(UP)^((At„)^)„ 

' 271^ ([/(t„) - /o]^)^ ((Ata)^), 


TT 

y 


(5.5) 


where in the last expression we have used the fact that 
since Op and Adg are small, |t/f — t/c| ^ Tmax- The two 


terms in (5.5) quantify the effects we predicted at the be¬ 


ginning of the section. The second term describes a loss of 
SNR due to the neutron star spin not being constant dur¬ 
ing the time spanned by a SFT pair, while the first term 
indicates a loss due to the mismatch between contribut¬ 
ing frequencies and the frequency of a single template. 
[In fact, the first term is just gf^f^ ([/(L) - /o]^)„-[ Note 
that we are free to choose the /g which maximizes the 
SNR for a given instantiation of spin wandering, which 
will be /o = (/(ta))^, so 

{[/(ia) - /o]")„ = ([/(L) - (/(L)).]^)„ (5.6) 

is the weighted variance of f{t) over the observing time. 

To get a quantitative estimate of the effects of spin 
wandering, consider a model where the neutron star spins 
up or down linearly with typical amplitude |/| drift, chang¬ 
ing on a time scale Tdrift where Tmax < Tdrift < ILbs- For 
simplicity, also neglect the impact of the weighting factor 

, so that (At^) « and (At^) « Then 




I/I 


drift 


(5.7) 


and 


([/(ia)-(/(io))«]^). 


< 


tct — Tmid 


Tdrift 

Fobs Tdrift i r|2 
~A I / I drift 


drift 


(5.8) 


Combining these results, we have 


So, in order to avoid a fractional loss in SNR of more 
than /i, one would need to limit the lag time to 


Tmax < min (l/ldriftVTobaTdrift) ' , ^^^l/ldrlft ^ 

(5.10) 

For example, if |/|drift = 10“^^ Hz/s, Tdrift = 10® s, 

Tobs = 1 Yr, and p = 0.1, the first limit is about 44,000 s 
and the second is 320,000 s. So in that case spin wander¬ 
ing would become an issue if T^ax ^ 12 hr. 

Note that this is somewhat less than the estimate 
AT ^ 3 day given in [6] . The source of this apparent 
discrepancy is a combination of the distinction between 
the coherent segment length AT and the maximum lag 
time Tmax, described in Sec. |IV B 3[ and the rough na¬ 
ture of some estimates used in |6]. That work compares 
the change in frequency |/|drift\/TobsTdrift/2 to the fre¬ 
quency resolution, which they give as ~ 1/AT. This is 
effectively an order of magnitude estimate, since it effec¬ 
tively assumes /i = 1, and also leaves out the numerical 
factor in l/y^g/y = •\/3/(7rAT). On the other hand, 
their frequency drift is the expected drift from the mid¬ 
dle of the run to the end; averaging the drift over the run 
gives an effective change of (| / 1 drift VTobsTdrift) /2 . In¬ 
cluding these three effects to do a calcul ation analogous 
to the one here would give a factor of tt a/5/3 « 4 reduc¬ 
tion on the estimated tolerable segment length to AT < 

2^3//^ (l/ldriftVTobsTdrift)”' « 62,000s « 17hr. Of 

course, the assumptions of |/|drift and Tdrift given above 
are uncertain and somewhat arbitrary, so our 12-hour 
number should also not be viewed as an exact constraint 
on the method. 


B. SFT length 


Most searches for continuous gravitational waves have 
used short Fourier transforms with a duration Tgft of 
30 min = 1800 s. The limiting factor which sets a maxi¬ 
mum on the reasonable Tgft is the accuracy of the linear 
phase approximation (2.15). 


If we consider higher order terms in the phase expan- 






















18 


sion, we have 


The effect of these corrections is to modify (2.21) to 


$(t(t)) « + 2tt fK{t - Ik) + - tK^ 

+ - txf + (tx)(t - txf + ... . 

(5.11) 




^tic+Tsft/2 


' —^sft/2 




^jtK) 

2 




F^A+-iF^A, 


ho(-l)'=e'^- + + 


dh(iK) / xt. 2 , . '$(iK) 


dt 


Io{KKk) + i ^ l2iKKk)Tsft +i J3(Kjc/c)yfft 

+ (5.12) 


where 


rl/2 

In{K) = / X'"e 
J-1/2 


”'■-‘2^'^^ dx = ( -^ ] 4^ sinc(K) 
' 27r 7 dK^ ^ ’ 


(5.13) 

Note that for even n, In{K) is real and even, while for 
odd n, it is imaginary and odd. 


We can then construct, as a replacement for (2.25) 


1 


^^K = y'' { — Io{l^Kk)hKk 

,—r 


fcG/Cjt 


hoe 


■^^^F^A+-iF^A^ Qk 


-if 


(5.14) 


Sk 


where 
Qk = 


^2 

“if 


.®(^if)v 
*— 


if02 T^ft + i 


3! ‘ 


rjio 

-'/COS sft 


$(iif) mx)? 


4! 


Sif04 


and 


SifOn = ^ IoiKXk)In{KXk) 
k^K,K 


(5.15) 


(5.16) 


As in Sec. |IIIB| we assume that the sum over pairs evenly 
and independently samples the fractional frequency offset 
kk from each SFT, which means we can replace Qk and 
Ql with 


(Qk).= (S^)+*^(So2)T;^ 


«>(iif) 

2 -1^02/ 4 sft 

^'{tK) [$(tif)]2' 

4! 8 


(5.18) 


{^04)T^tt 


where the fact that Is^k) is odd in k means that the 
average (Soa) vanishes. 

Now, 


Re 


(QifQLrifi) = Re (QxQl) ReFft-i—Im (QxQl) IniF/f-L 


Re (QxQl) ^ iQxQl) 


KL 

(5.19) 


We assume that the impact of the second piece is smalj^ 
and focus only on Re (QxQl), which leads to a fractional 
loss of SNR of 


1 - 


E[p] 


/■='2\2 


)^-{ReiQKQl)) 


_ ( i^x) + i^l) (S04) (^if^L) (So2)2 


/" 2 \ 


rpA 

^2\2 I ^sft 


The expectation value (3.7) of the statistic thus becomes, 
including the correction for higher phase derivatives and 
finite SFT length, 

E[p\^Nhl2 Y. fTL^e(QKQlfKL) (5.17) 

KL^V 


(5.20) 


In particular, it is suppressed by averaging non-positive-definite 
antenna patterns, although the same combination is the source 
of systematic errors in parameter estimation. 
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Differentiating (4.15) gives 


= 27r/od_R- — 


PL 


foQp sin 


27r 

pT 


orb 


i^K ^asc) 


(5.21) 


We can neglect the first term, since the acceleration due 
to the Earth’s orbit is 0(10“^^ s“^) and that due to the 
Earth’s rotation is 0(10“^° s“^). In comparison, for Sco 
X-1, 


27r 


orb 


= 1.23 X 10"® s"^ 


(5.22) 


If we assume, as in the metric calculation, that the aver¬ 
age over pairs evenly samples the orbital phase, then 




p4 

-"orb 


Using the identity 


sin A sin i? = ^ [cos(A — B) — cos(A -I- B)] 


(5.23) 


(5.24) 


we can calculate 
27r 


sm 


-^orb 


(tiC ^asc) 


sm 


27rAtn 


COS ■ 


— ( COS 


27r ^ \ 

iJ'L ^asc) 

-torb 

tasc) 


-forb 

SO the fractional loss in SNR is 

E{P] 

P M ideal 

( (S 04 ) (So2)^ 


-forb 


(5.25) 


1 - 


p4 

■^orb 




COS • 


27rAio 

P. 


orb 


rjiA: 
^ sft 


(5.26) 


The f actors (E 04 ) and (E 02 ) can be calculated by using 
(5.16) along with 


h{K) = 


SlUTTK COSTT/s: 


SlUTTK 


Airti 2{7 tk)‘^ 2{7tk)^ 


(5.27) 


and 


hin) = 


sin TTK cos TTK 3 sin ttk 3 cos ttk 3 sin ttk 


IhTTK 4(7rK)2 4(7r/c)3 2(7rK)^ 2(7rK)® 

(5.28) 

and averaging numerically over k given the number of 
frequency bins included. In Table |IV[ we show the two 
coefficients appearing in (5.26), for various choices of the 
number m of included frequency bins (see also Table 0 - 


Note that for the cross-correlation search, choosing 
shorter SETs does not directly impact the sensitivity. Eor 
the same allowed lag time, searches with different SFT 
lengths should have approximately the same sensitivity. 


TABLE IV. The c oefficients (Eo4)/{h 2) and (Eo2)V(2^>^ 
appearing in (5.261, for various choices of the number m 
of included frequency bins, where (Son) is the mean value 
of Eo„(k) = Y}LLP-4)/•2^ P s)Irt{n + s), averaged 
over — I < AC < and In{K.) is defined in (5.131 with 
/o(At) = sine AC = B^k) are given by 

(5.271 and (5.28). Note that the value of (H^) = {loo) is 
tabulated in Table El 


m 

1 

2 

3 

4 

5 

6 

(So4)/{h2) 

0.0107 

0.0056 

0.0086 

0.0042 

0.0099 

0.0052 

0.0100 

0.0055 

0.0106 

0.0059 

0.0108 

0.0060 


We can see this by considering the SNR for a given signal 
amplitude /iq, for example, from (3.16). Since 


■nave -p 
^ KL = ^ KL 


2'rsft 


(5.29) 


the quantity (r^£)^ inside the sum is proportional to 
(Tsft)^. However, for a fixed maximum time lag Tmax, 
the number of terms in the sum will be proportional to 
(Tsft)"^ and the resulting expected SNR will be approx¬ 
imately independent of Tgn- (For example, halving the 
SFT length will mean each SFT pair contributed one- 
fourth as much to the sensitivity, but will double the 
number of SFTs and thus quadruple the number of SFT 
pairs.) 

On the other hand, by increasing the number of SFT 
pairs, using a shorter SFT length will mean increasing 
computing cost at the same Tmax- If the computing 
budget is f ixed, the sensitivity gained by reducing the 
mismatch (5.26) will be offset by the loss of sensitiv¬ 
ity, in the form of a lower resulting from a 

smaller Tmax- Following the reasoning in Sec. IHT^ if 
the computing cost scales like the number of templates 
(which scales like Tj^ax) times the number of SFT pairs 


(which scales like T)nax7obs7)jft ); then the overall sensi¬ 
tivity for a fixed observing time Fobs scales like ^ 

and therefore the restriction at constant computing cost 
2 

will be Tniax cx . Since the sensitivity scales with 
the square root of the number of SFT pairs, we have 

E cx and 




where 


. 8^6a2/(Eo4) (Eo2)V 27TAE 

~ —zza I , ^ \ COS 


PL V(s^) 


/'='2\2 


Eorh 


(5.30) 


(5.31) 


is the mismatch scaling appearing in (5.26) ^ The sen- 


21 Of course A still depends on Tnax through ^cos ^; but 
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Frequency (Hz) 


FIG. 4. The optimal SFT length Taft, defined in ( 5.34| and 
(5.311, as a function of frequency, for a signal with the most 
likely orbital parameters for Sco X-1, as given in Table |I] as¬ 
suming that d = 3, i.e., the density of points in parameter 
space grows as the third power of the coherence time Tmax. 
This is appropriate for a search over, e.g., freqnency fo, pro¬ 
jected semimajor axis Up, and time of ascension tasc (when the 
uncertainty in the period Porb is small enough that a single 
value may be assumed), in the case where Tmax is small com¬ 
pared with Porb. The solid line represents a more optimistic 
scenario where the average cosine appearing in the second 
term of ( |5.31[ ) is approximately unity, which should also be 
the case if Tmax <C Porb. The dashed line represents a worst- 
case scenario where the average is approximately zero. The 
optimal SFT length maximizes the expected SNR in (5.301 
and represents a balance between two competing effects: if 
Tsft is too large, phase acceleration will lead to a loss in SNR 
compared to the ideal formula (3.191; if Taft is too small, the 
large number of SFT pairs in the computation will lead to a 
restriction on the possible Tmax achievable at fixed computing 
cost, and reduce the ideal SNR itself. 


sitivity at fixed computing cost is thus maximized when 
1 _ (4d + 5)AfX\ = 0 (5.32) 

i.e., when the mismatch due to SFT length is 

fi = (5.33) 

The corresponding optimal SFT length is 

Taft = ([4d+ 5]A)-i/Vo”'^" (5-34) 

For example, if d = 3, ^ « 0.059. In Fig. H we 

show this optimal SFT length for d = 3, using Up = T44 s 
and Porb = 68023.70 s (the most likely values for Sco X- 
1). The solid line shows the most optimistic scenario, 

in which « 1 (which will be the case for 

Tmax Porb) and the dashed line shows the most pes¬ 
simistic scenario, in which the average goes to zero. 


if Tinax is small compared to Porb? which we are assuming in 
the scaling of number of templates with Tmax, this average is 
approximately unity. 


VI. CONCLUSIONS AND OUTLOOK 


In this paper we have explored details of the model- 
based cross-correlation search for periodic gravitational 
waves, focusing on its application to signals from neu¬ 
tron stars in binary systems (LMXBs) and Scorpius X-I 
in particular. We have carefully considered the impact 
of spectral leakage (in Sec. Ill B) and the implications 
of unknown amplitude parameters (in Sec. Ill CI on the 
sensitivity of the method. We have also produced ex¬ 
pressions for the parameter space metric of the search 
(in Sec. IVB I, at varying levels of approximation, and a 
systematic offset in the parameters of a detected signal 
related to the unmeasured inclination angle of the neu¬ 


tron star to the line of sight (in Sec. IVA). In Sec. VA 


we estimated the effects of “spin wandering” caused by 
deviations from equilibrium in the torque balance config¬ 
uration, and in (V B) we consider the appropriate SFT 
duration needed to avoid significant loss of SNR due to 
unmodeled phase acceleration. 

We have shown (in Sec. HID I that the method pro¬ 
duces an improvement in strain sensitivity over the di¬ 
rected stochastic search method which inspired it; this 
is roughly proportional to the fourth root of the product 
of the coherence time of the model-based search and the 
frequency bin size for the stochastic search. A mock data 
challenge [32] has been carried out by comparing the per¬ 
formance of the available search methods, including the 
model-based cross-correlation search, on simulated sig¬ 
nals injected into Gaussian noise. As reported elsewhere 
PSI 132] , the cross-correlation search is the most sensitive 
one currently implemented. 

To give an estimate of expected sensitivity for data 
from detectors such as Advanced LIGO and Advanced 
Virgo, it is necessary to make some suppositions about 
the parameters of the search, especially the time T^ax 
over which SFTs are coherently cross correlated. Since 
this drives both the sensitivity and computing cost, the 
choice of Tmax will depend on available computing re¬ 
sources, and will likely vary with frequency in order to 
optimize the distribution of computing resources where 
they can be most effective. In [55], we performed searches 
with 9 min < Tmax ^ 90 min for a range of frequency 
bands covering a total of 500 Hz distributed in /g G 
[50,1455] Hz, using moderate computational resources. 
On the other hand, in Sec. |V A[ we considered spin wan¬ 
dering effects which might lead to a significant loss of 
SNR for a search with Tmax ^ 12 hr for a one-year obser¬ 
vation. 

In Fig. we show the projected sensitivity (3.25) of 
a search using one year of data, either from the two ad¬ 
vanced LIGO detectors in Hanford, WA and Livingston, 
LA, or from the two advanced LIGO detectors plus the 
Virgo detector in Cascina, Italy, all operating at their 
projected design sensitivity. We show the sensitivity of 
three hypothetical searches, with Tmax = 6 min, 60 min or 
600 min = 10 hr, and compare the observable Hq (at a 5% 
false dismissal probability, assuming a single-template 
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FIG. 5. Expected sensitivity (3.251 for a search of one year 
of coincident data from either the two LIGO detectors (la¬ 
beled HL) or the three LIGO-I-Virgo detectors (labeled HLV), 
at design sensitivity. The value plotted is the observable 
ho at 5% false dismissal probability, assuming an overall 
false alarm probability of 5% and a trials factor of 10® for 
a single-template false alarm probability of 5 x 10“^° (i.e., 
see Sec. IIIG and Table Till. The three curves in each set, 
are, from top to bottom, for Tmax = 6 min, 60 m in and 10 hr. 
They are compared to the signal strength (6.21 predicted by 
the torque balance argument |12|. 


false alarm probability of 5 x 10 corresponding to an 


overall 5% false alarm proba bility and a t rails factor of 
10®, as described in Sec.llII Gland Tablelllll). For compar¬ 


ison, we show a representative signal strength predicted 
by the torque balance argument [HI [H]. By assuming 
that the spin-down torque due to gravitational waves is 
balanced by the spin-up torque due to accretion, esti¬ 
mated using the observed x-ray flux, it is possible to es¬ 
timate the strength of the gravitational-wave signal as a 
function of the neutron star spin frequency [H]: 


/ig 


3 X 10"^^ 

X ( ^ 

^ I 10km 



Fx 

erg cm“2 

M \ 
IAMq J 



1/2 


12s \-V2 

300 Hz/ 


( 6 . 1 ) 


The spin frequency of Sco X-1 is unknown, but Vg val¬ 
ues inferred for other LMXBs from pulsations or burst 
oscillations range from 50 Hz to 600 Hz, so we consider 
the sensitivity over a wide range of GW frequencies. 
For Sco X-1, using the observed x-ray flux Fx = 3.9 x 
10“^ ergcm“^ s“^ from [H], and assuming that the GW 
frequency /o is twice the spin frequency Vg (as would be 
the case for GWs generated by anisotropies in the neu¬ 
tron star), the torque balance value is 


ho 


3.4 X 10"^® 


Vs 

300 Hz 



( 6 . 2 ) 


which is the reference curve plotted in Fig. We see 
that for a three-detector, one-year analysis, a signal at 


FIG. 6. Relative scaling of expected computing cost per log¬ 
arithmic frequency interval for a search of one year of coinci¬ 
dent data from either the two LIGO detectors (labeled HL) 
or the three LIGO-|-Virgo detectors (labeled HLV). The three 
curves in each set, are, from top to bottom, for Tmax = 10 hr, 
60 min and 6 min. The calculation assumes that the comput¬ 
ing cost scales with the number of SFT pairs times the number 
of points in parameter space. It also assumes that the optimal 
SFT length Tsft given by (5.341 and (5.311 has been chosen 
at each frequency, and that we are searching over frequency 
and two orbital parameters. The approximate scaling is as 
/o/max, so for instance a Tmax = 60 min search from 100 to 
200 Hz would consume the same resources as a Tmax = 6 min 
search from 1000 to 2000 Hz. For reference, the mock data 
analysis in m, which was accomplished in approximately 
20,000 GPU-days, covered a set of roughly logarithmically- 
spaced frequency bands totaling 250 Hz spread from 50 Hz to 
1375 Hz at a range of Tmax values from 9 to 90 min. 


the torque balance limit should be detectable for 30 Hz ^ 
/o 300 Hz with T^ax = 60 min (which is already com¬ 
putationally manageable at most frequencies), and if one 
could increase to T^ax = 600 min through algorithmic 
improvements, programming optimization, and/or appli¬ 
cation of additional resources, that range could be broad¬ 
ened to 20 Hz ^ /o ^ 500 Hz. The best-case ho sensitiv¬ 
ity of 5 X 10“^® for the 60 min search is consistent with 
the results of the Sco X-1 MDC [H], where a cross¬ 
correlation search with 9 min < T^ax < 90 min was able 
to detect signals with hg ^ 5 x 10“^®. 

The choice of T^ax will in part be constrained by com¬ 
puting cost; in Fig. we show the approximate rela¬ 
tive computing cost scaling for the six searches consid¬ 
ered (one year of data from either the two LIGO de¬ 
tectors or the two LIGO detectors and Virgo, with a 
maximum allowed lag time of Tmax = 6 min, 60 min or 
600 min = 10 hr. The computing cost is assumed to be 
proportional to the number of SFT pairs times the num¬ 
ber of parameter space points to be searched, and we 
plot the relative cost per logarithmic frequency interval. 
We also assume that at each frequency the SF T leng th is 
chose n to be the optimal SFT length given by (5.34) and 
(5.31). Roughly speaking, the number of SFT pairs will 
scale as /gTmax (since the optimal SFT length scales as 




































































22 


_ /2 

'Tmax ), and the density of templates in parameter space 
will scale as /oT^ax- The density of points per logarith¬ 
mic frequency interval introduces another factor of /o, 
so the quantity plotted, cost per unit frequency interval, 
scales approximately as /oT,^ax- This means that, for ex¬ 
ample, a Tniax = 60 min search from 100 to 200 Hz would 
consume the same resources as a T„iax = 6 min search 
from 1000 to 2000 Hz or a Tmax = 600 min search from 
10 to 20 Hz. 


Finally, we consider one possible avenue for enhance¬ 
ment of the cross-correlation method. As explained in 
Sec. HI A the fact that we filter with means that 
the method provides an estimate of /iq®, a function of ho 
and cos 6 dehned in (3.10), rather than Hq. If we had a 
method of independently estimating cos t, or in fact any 
other combination of ho and cos t besides /iq®, we could 
obtain a better measurement of ho- In EB, a method 
was proposed to obtain estimates of hoA+ and hoAx , 
but a more effective procedure would se em to be adding 
a second statistic which uses ir™£ [see (2.32b)] in place 
of and therefore observes the quantity h'oA+Ax ; be¬ 
tween this and the original /i^® estimate, we would be able 
to disentangle ho and cos i. This prospect bears further 
investigation. 
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51,Sh - Ik) 

+Tsft/2 


JtK — Tsft/2 \ ^Sft / 
. 1/2 


-i27T{fk-fK)it-tK) 


= rsft / = (A2) 

Jiji 

with KKk = ifk — /if)Tsft as befor e. The noise contribu¬ 
tion is modified by replacing (2.8) with 

A[nMi1«^ifL7riTsft^ (A3) 


where 


Appendix A: Effects of Nontrivial Windowing 


1. General formulation 


As noted in Sec. H A the construction of Eourier trans¬ 


formed data is often done with a window function w{9), 
as in (2.3), as opposed to the unwindowed (or nearly- 


rectangularly-windowed) data considered in the main 
body of the text. This appendix considers the impact 
on the search method and its sensitivity of using a non¬ 
trivial window function, which is investigated in greater 
detail in m- 

The use of windowing for Fourier transforms affects 
the expected signal and noise contributions to the data. 


Ike = 


i-ir 


k—e r°° 


' 6lJfk-f)5£^ifi-f)df 

-‘-Sit J —oo / A A\ 

1-1/2 (A4) 

= (-l)'=-M e*2-('=-^)«[u;(0)]2d0 

J-1/2 

Note that the diagonal elements of this matrix are equal 
to the mean square of the window function: 


/•1/2 _ 

lkk= [wie)fd9^w^ (A5) 

■ 1 - 1/2 


If we define 


^Kk — -^Kk 




(A6) 
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as in (2.91, we have 


E [^'ku] — Mk/c 


7K 


?K 


and 


E [{zlu - ^^lk){zLl - = Skl 7S • (A8) 



1.0 


0.8 

2'rsft 


Sk 

^ 0.6 

(A7) 

T 


0.4 


0.2 

(A8) 

0.0 


We then modify (2.23) to 

E (A 9 ) 

keKK k'&KK 

where {( 7 “)^/} are the elements of the matrix inverse of 
{7feJ, and 



FIG. 8. Specific versions of the general Tukey window wp(6) 
as defined in ( |A12[ ): the rectangular window = wo{d), 

a canonical (/3 = |) Tukey window tCi/ 2 (^), and the Hann 
window = wi{6). 


'^W 


Y fcG/Cjf fc'GKjc 

(AlO) 


ensures that the normalization (2.24) holds as before. 


Then the derivation proceeds as before, with S^_replac- 
ing Ek, and, in particular, the expected SNR (3.16) be¬ 
comes 


E[p]^{hfr{{^nYj2 E (r??i) (All) 

V KLeV 


2. Results for specific windows 

We now c onsider the consequences of the modifi¬ 
cation (All) by investigating the form of ^^{k) = 


'Esh^T^tS'^/'Esh) defined in (A2) and 7 ^^ defined in (A4) 


for specific nonrectangular window choices. We consider 
the general family of Tukey windows, defined using an 
adjustable parameter 0 < /3 < 1 by 




(^1 - cos I (20 -I-1)^ 

wpiO) = <( 1 - ( 1 ^) <0<(- 

(1 - cos f (20-1)) (1/)<0<1 

(A12) 

The general form of the Tukey window is illustrated in 
Fig.0 This family includes at its extremes the rectan¬ 
gular window (/3 = 0) and the Hann window (/3 = 1). In 
practical applications it is also common to use a Tukey 
window with a small finite parameter /? <C 1 rather than 
a pure rectangular window. These two specific cases are 
shown in Fig. 8 along with a Tukey window with /3 = 1. 


into (A2) to obtain 


CY(«^) = EincK-h ^(1 -/3)sinc(K[l -/3]) 


/3 . , 

sm I TTK 


1-i 


1 — /3^\ . f 1 13 k, 

sine 1 —-— — sine 


(A13) 


the “interesting” values of /3 also have somewhat sim¬ 
pler explicit forms. For the rectangular window (/3 = 0), 
which was considered in the main body of the paper, we 
have 


^“(k) =f®“‘(K) =sincAC ; 


for the Hann window (/3 = 1), we have 


(A14) 


We can insert the general form of wp{9) from (A12) 


^“(k) = i sine K-l-i sinc(l— k)-|-^ sinc(l-l-K) ; 

(A15) 

and for the canonical (/3 = ^) Tukey window, we have 

er/2(«) = e^"'^"(«) 

= - sine K -sinc(2 -f k) -sinc(2 — k) 

2 4 ' ^ 4 ' ' 

-b i sine Y + l sine (l + |) + ^ sine (l - ^) (A16) 


We plot these three functions in Fig.[^ 

To evaluate the factor of ((S™)^) appear ing in (All), 
we need to construct the matrix { 7 ^^} via (A4). Substi- 
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FIG. 9. The window leakage function = 

(^^/Tsft) defined in (A2l for the windows shown in 
Fig. i The explicit formulas are given in ( |A14[ ) for the rect¬ 
angular window, (A16l for the canonical (/3 = |) Tukey 
window, and (A15l for the Hann window. Note that the 
version for rectangular-windowed data is just = 

(fv/Tsft) = sinc(K) = > which is the finite-time 


delta function plotted in Fig. 



FIG. 10. The leakage factor appearing in (All I for 

a search using between one and six bins from each SFT, as¬ 
suming a general Tukey window from the family (A12l. We 
see that, for any number of bins, the most sensitive search is 
when /? = 0, i.e., for rectangular windows. In particular, when 
a single bin is used from each SFT, we have (H^) = 0.774 for 
rectangular windowing (/3 = 0), ((H§’)^) = 0.699 for a canon¬ 
ical (/3 = I) Tukey window, and = 0.601 for Hann 

windowing (/3 = 1). Note that the /3 = 0 value on each curve 
is just the corresponding “cumulative” number from Table [III 


tuting (A12) into (A4), we can find 


{rp)M = - P) Sine [{k - £)(1 - /3)] 

+ ^/3sinc[(fc-7)/3] 

— -/3 sine [(fc — i)P — 1] — -/? sine [(fc — £)f3 + 1] 

+ sine [{k - i)(3 - 2] + sine [{k - i)l3 + 2] . 

Id Id 

(A17) 

We can see that, for the rectangular case /3 = 0, we get 
[lo)kt = Ske as before, while for the Hann case /3 = 1 we 
have 




The diagonal elements for general /3 are 
il^)kk = 1 - ^/? = 


1 

16 


^k,e+2 ■ 
(A18) 


(A19) 


as in (A5). This means that, in the special case where 


the set of bins ICk f rom each SFT is just the “best bin” 
kx defined in (2.18), the matrix just has a single 


element r = 1 — |/3, and 

'fcjffcic 


= 




1 - 


(A20) 


where ^^(k) is defined in ( |A13| ). In general, though, we 
need to invert the matrix ( A17| ) and then average 
defined in (AlO) over possible values of kk- We plot the 
results in Fig. m as a function of /3, for cases where we 
take the “best” m bins from each SFT. We see that, for 
any number of bins, ((S“')^) is a maximum for /3 = 0, i.e., 
rectangular windowing. The (3 = 0 values are just the 
“cumulative” entries from Table |ll] for the corresponding 
number of bins. Specifically, for the single-bin case, when 
(3 = 0, we have (S^) = 0.774 (as seen in the m = 1 entry 
of Table [n| , when (3 = ^, we have ((2'^'^^®^)^) = 0.699, 
and when (3 = 1, we have ((2^“®)^) = 0.601. These 
values also appear in m, which explains in more detail 
the relevant phenomenon. While the dropoff from the 
maximum value of to its average value is great¬ 

est for rectangular windowing, the maximum value and 
the average value are also greatest for the rectangular 
window. 

A common approach to handle the loss of signal as¬ 
sociated with Hann-windowed data is to divide the data 
into overlapping Hann-windowed data segments, as in 
[l8] . For the present search, however, it is easier just 
to include more bins from the rectangularly windowed 
Fourier transform, if desired, to increase the sensitivity 
of the search. The only drawback to that is a slight in¬ 
crease in computational time, but this increase is much 
smaller than what would arise from almost doubling the 
number of SFTs by the use of overlapping windows. 
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FIG. 11. Eigenvalues {lok} of the weights matrix W defined in(2.35l for two scenarios. On the left, we show one day of 
observation with the LLO, LHO and Virgo detectors, assuming equal sensitivity, with Taft = 900 s and Tmax = 3600 s. On the 
right, we show one year of observation under the same conditions, constructed as the union of 365 such days, spread throughout 
the year. In both cases, the start and end of each day include data gaps of 900-1800 s, randomly and independently generated 
for each detector. 


Appendix B: Probability distribution for 
cross-correlation statistic in Gaussian noise 


In this appendix, we consider the detailed statistical 
properties of the cross-correlation statistic (2.36) in the 


presence of Gaussian noise. If the noise contribution to 


XKk is Gaussian, the definitions (2.9) and (2.23) imply 


that z—/X is a circularly symmetric Gaussian random vec¬ 
tor |35j with zero mean, unit covariance and zero pseu¬ 
docovariance, as described in (2.26). If {ujk} and {vx} 
are the eigenvalues and eigenvectors, resp ective ly, of the 
Hermitian weighting matrix W defined in (2.35) , so that 


W = ^ 


(BI) 


K 


then the statistic is 


K K 


V^Z 


(B2) 


The conditions Tr(W) = 0 and Tr(W^) = I imply that 
— 0 and = example of 

the typical form of the eigenvalues, we present in Fig. |11| 
two typical sets of eigenvalues, one assuming a day-long 
observation with three detectors, assuming Tgh = 900 s 
and Tniax = 3600 s, the other combining 365 such obser¬ 
vations with randomly staggered starting times to simu¬ 
late a year-long observation, assuming LIGO Livingston, 
Hanford and Vir go detectors with identical and station¬ 
ary noise spectrap^ 


Note that since , a matrix made of the 

{r^£} has the same eigenvalues as one made of the If 

the noise PSDs are (approximately) the same for all SFTs, it is 
also equivalent to using the eigenvalues of a metric made of the 


Each v^z is an independent circularly symmetric 
Gaussian random variable with zero mean and unit vari¬ 
ance, which means its real and imaginary parts are inde¬ 
pendent Gaussian random variables with mean zero and 

is ^ times a X^(2) random vari¬ 


variance Thus 


vJfZ 


able, i.e., it is an exponential random variable with unit 
rate parameter. The characteristic function is thus 


(pK{t) = E 


vL 


1 — it 


(B3) 


which means that the characteristic function of the cross¬ 
correlation statistic is 


= E 


exp 




K 


1 

riif (1 - 


(B4) 


This allows a straightforward computation of the exact 
probability density function for the statistic p as 


f{p\K = 0) 


E 




K,ujk>0 


K,ujk<0 


p > 0 
p < 0 


(B5) 

which is a a mixture of exponential distributions. To get 
the false alarm probability a at a threshold p*** > 0, we 
calculate 


a = P{p > p“'|ho = 0) 


/(p|/io = 0)dp 


= E 

u;jc>0 




(B6) 


The problem with this expression is that the denominator 
can get very small, and the signs of the terms alternate. 
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FIG. 12. False alarm probabilities for the cross-correlation 
statistic in the day-long and year-long scenarios considered in 
Fig. 11 using the explicit formula (|B9[) as well as numerical in¬ 


tegration of (BIO I, along with the probabilities we would get if 


we assumed the statistic to be Gaussian. For a day-long obser¬ 
vation (with three detectors, = 900 s and T^ax = 3600 s), 
both methods give comparable results, but the Gaussian ap¬ 
proximation is invalid for single-template false alarm proba¬ 
bilities below about 10“^. Note that for large signal values, a 
single exponential term dominates. For a year-long observa¬ 
tion, practical calculation with ( |B9[ ) is impossible due to un¬ 
derflow issues. The numerical integration of (BIOI becomes 
unstable for false alarm probabilities below 10“^^, but not be¬ 
fore quantifying deviations from the Gaussian approximation 
even for a year-long observation. 


To see this, assume that we have ordered the eigenvalues 
so that 


UJn > UJN-1 > ■ ■ ■ > UJKo > 0 > OJKo-I > ■ ■ ■ > UJi (B7) 


Then 


no 


L^K 
= (-1)^-^ 


UJk 


'K-1 


nd 
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LOk 
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n (1 


UJk 


‘K-1 


nd 
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LOk 
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(B8) 


and the false alarm probability is 


N 

a= ^ (-1) 

K^Ko 

rK-1 


N-K^-P^^/uk 


X 


nXS) 


-1 r 


N 

n 

Lz,=if+i 


UJk 


- 1 


-1 


(B9) 


The last two factors can be very large, and are larger 
when the eigenvalues are closer together. (Recall that N 
is the number of SFTs, which is approximately TLbs/Tsft, 
so there are many factors appearing in the product.) 

Given the numerical problems with the exact false 
alarm probability (B9) when the number of SFTs is large. 


it is sometimes necessary to use an alternate approach. 
We can perform a calculation analogous to that in [18], 
based on the method of ISSIIST]. This uses the Gil-Pelaez 
expression |38j to construct a cumulative distribution di¬ 


rectly from the characteristic function (B4) according to 


1 1 
“=2 + ^ 


Im 


= -itp 


dt 

~T 


(BIO) 


We can then find the false alarm probability by numerical 
integration of (BIO). Results of both of these methods 
are shown in Fig. m for the two scenarios considered in 
Fig.[TT] Both methods produce consistent results for a 
day-long observation, and illustrate deviation of the false 
alarm probability from the Gaussian value for > 2. 


For the year-long observation, explicit evaluation of (B9) 


is impossible because of underflow in the cancellations, 
but numerical integration of (BIOI works until the false 
alarm probability goes below 10“^^ or so. False alarm 
probabilities are considered in detail for a wider range of 
observing scenarios in |55|. 
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